QQ-onia package: a numerical solution to the Schrodinger radial equation 
for heavy quarkonium 

Juan-Luis Domenech-Garret * ^ 
Miguel- Angel Sanchis-Lozano ^ ^ 

'^Departamento MACS, Ffsica Aplicada. Universitat de Lleida. 
Av. Alcalde Rovira Roure 191, Bldg-3, E-25198, Lleida (Spain) 

'^Institute de Fisica Corpuscular (IFIC) and Departamento de Ffsica Teorica, 

Centre Mixto Universitat de Valencia-CSIC 

Dr. Moliner 50, E-46100 Burjassot, Valencia (Spain) 

This paper presents the basics of the QQ-onia package, a software based upon the Nu- 
merov 0{h^) method which can be used to solve the Schrodinger radial equation using a 
suitable potential V{r) for the heavy quarkonium system. This package also allows the 
analysis of relevant properties of those resonances such as the square of the wave functions 
at the origin, their corresponding derivatives for / 7^ states, typical heavy-quark veloc- 
ities, and mean square radii. Besides, it includes a tool to analyze the spin dependent 
contributions to the heavy quarkonia spectrum, providing the splitting of n^Si —u^Sq, as 
well as the n^Pj — n^Pi energy levels. Finally a simple software implemented in QQ-onia 
to compute El transition rates is presented. 

PROGRAM SUMMARY 

Program title: QQ-ONIA PACKAGeI. 

Manuscript title: QQ-onia package: a numerical solution to the Schrodinger radial equa- 
tion for heavy quarkonium. 

Manuscript Authors: Juan-Luis Domenech-Garret; Miguel-Angel Sanchis-Lozano. 

Programming language: PAW (Physics Analysis Workstation). 

Operating system(s) for which the program has been designed: Windows-XX and Unix 
(Linux). 

Keywords: Heavy Quarkonium potential; Wave function at the origin. 
PACS: 14.40.-n; 12.39.-x; 14.40.Gx 



* Corresponding author E-mail: domenech@macs.udl.es 
'I' E-mail:Miguel. Angel. Sanchis@uv.es 

•^Program Author and Copyright: Juan-Luis Domenech-Garret 



2 



Nature of the problem: Software to solve the Schrodinger radial equation using a suitable 
potential V{r) for the heavy quarkonium system, allowing to perform spectroscopy. 
It also allows the analysis of relevant quantities of those resonances such as the 
square of the wave functions at the origin, their corresponding derivatives for Z 
states, typical heavy-quark velocities, and mean square radii. The package is a (user- 
friendly) multipurpose tool for dealing with different heavy quarkonium systems, 
providing a way to study the influence of a given potential on a scries of relevant 
physical quantities, by either varying parameterized values of a well-known potential 
form, or by including new terms. 

Solution method: Based upon the Numerov 0{h^) Method, we perform a matching proce- 
dure to the reduced wave function at the cut point. We also perform a normalization 
technique for for these wave functions taking into account the different domains when 
we use a Numerov backward-forward technique. In the case of / > 2 we present a 
way to find the corresponding derivatives at the origin by only calculating the re- 
duced radial wave function and first derivative. When estimating the heavy quark 
velocity, we introduce an additional way to compute this quantity from the virial 
theorem. The calculated reduced wave functions and radial wave functions at the 
origin are later used to obtain the heavy quarkonia nL splitting and El transition 
rates. 

Additional comments: Using Windows, to optimize the edition of the files, please, open 
it with MFC-WORDPAD. 

Running time: It depends on the choice of the r range, and the number of energy steps. 
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1. Introduction 

Since the discovery of the charmonium and bottomonium families, much efforts have 
been spent over 30 years trying to understand the nature of heavy quarkonium, (as a 
summary see for example pQ); in the meantime various numerical tools have been created 
in an almost "ad hoc" fashion with the aim of extracting important information about 
their properties, such as their masses and partial widths. 

This package was developed with the aim of providing a multipurpose (user-friendly) 
tool for dealing with different heavy quarkonium systems, providing a way to study the 
influence of a given potential on a series of relevant physical quantities, by either varying 
parameterized values of a well-known potential form, or by including new terms. 

The QQ-onia package handles the heavy quarkonium system within a non relativistic 
framework, solving the Schrodinger radial equation [SRE) with an appropriate potential 
for heavy quarks/antiquarks. The basic reason for this choice derives from the Quark 
Potential Model [2j, which establishes a low value for the expected square velocity of the 
quark for these heavy resonances (y"^ ~ 0.1 for bottomonia and ~ 0.3 for charmonia). 
Besides, there is also another reason from a dual ultra-relativistic picture to provide a non 
relativistic treatment for heavy quarkonium [5]. These low velocities << 1 were also 
responsible for the success of the Non Relativistic QCD (NRQCD) ([1], [S]), a rigorous 
effective theory for strong interactions deriving from first principles. 

The paper is organized as follows: the first section explains the basics of the heavy 
quarkonium Non Relativistic potential and its spin-depedent part, according with the 
Breit-Fermi Hamiltonian. Later we explain the underlying foundations to our way of 
solving the Schrodinger equation. Although the arguments presented are already well- 
known, we consider this introduction necessary to facilitate the understanding of the rest 
of the article. We then look into the specific details about the calculation of the wave 
functions of heavy quarkonia in our code. 

We initially focus on a procedure for matching the reduced wave function at the cut 
point and its normalization when we use a Numerov backward-forward technique. We 
also explain how to extract useful quantities in order to analyze the heavy quarkonium 
system such as mean square radius, and heavy quark velocity (introducing an additional 
way to calculate this quantity from the virial theorem). For / = states we also focus 
on square of the wave functions at the origin, and, for I ^ states, we examine their 
corresponding derivatives; in the case of Z > 2 we present a way to find these values by 
only calculating their corresponding reduced radial wave function and first derivative. It 
follows an explanation on the spin-dependent terms in the bottomonia case and, finally, 
we focus on the calculation of the El transition rates. 

The software will be explained in more detail in the second part of the paper, where 
we also present several results obtained using the programme to illustrate the procedure 
and foreseen precision. 
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Several potentials for the heavy quarkonium system are given as examples. Our first 
(main) choice is the well-known funnel or Cornell potential, i.e. a Coulomb plus Lineal 
static potential (CpL) [6j, 

Cpash 

V [r) = ar 

r 

with string tension cr and strong coupling constant Og] Cp is a colour factor. 

Nevertheless, other static and non-static potentials will also be used in this paper in 
order to illustrate some results. Note that we will focus on the equal masses. Thus one 
can write for bottomonium rriQ = mg = = m with reduced mass n = wiq/2. 

2. Physical Bases of the package 

2.1. Heavy quarkonium potential and spin-dependent terms 

If VjsfR stands for the Non Relativistic potential, one can split it in two terms consisting 
of a vector (V) and a scalar (S) contribution [7f8] 

VNR{r) = Vv{r) + Vs{r) (1) 

In our example with a funnel type potential (being k = Cpagfi) 

Vv{r) = —k/r ; Vs{r) = ar 

In accordance with literature (e.g. [7J), additional terms have been included in the 
potential Ynr to take into the account the spin orbital and the spin-spin interactions, 
causing the splitting of the different mass levels. The additional potential reads [B| 

V{r) spin- dependent = VlS + KsS + (2) 

where Vis, Vss, and Vp are the spin-orbit, the spin-spin, and the tensor terms, respec- 
tively. The spin-orbit term, in the equal quark masses case, is 



VLs{r) 



3 -VV(r) - ^Vsir) 



2 r I dr dr 



(3) 



where L is the relative angular momentum of the constituents (1 and 2), and S is the total 
spin of the bound state, S = Si -|- S2 ( with J = L -|- S ); ((L • S)) for different j and / 
values is shown in Table 1. 



Table 1 

((L • S)) coefficients ((L ■ S) = if / = or 5 = 0). 



j value 


+ 


/ 


(/-I) 


(L-S) 


/ 


-1 
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The spin-spin term can be written as 
2 (Si ■ S2) 



Vssir) 



A( Vv(r) 



3 r 

where the (Si ■ S2) coefficients take the values: 

(—3/4) for the spin-singlet case {S = 0), and 

The tensor term can be written as 
1 r/ 

r ar dr"^ 



VAr] 



12 rri- 



is 



12 



(4) 



-1/4) for the spin-triplet case {S = 1). 



(5) 



where (5*12) is the spin-dependent factor (for / 7^ and S = 1), shown in Table 2 for 
different j and / values. 



Table 2 



Spin-dependent Su factor. ( (5'i2) = 0if/ = 0orS' = 0) 



j value 


(^ + 1) 


/ 


(/-I) 


(^12) 


- 2 //(2/ + 3) 


2 


- 2 (/ + l)/(2/- 1) 



2.2. Schrodinger radial equation 



Basically, our code has to solve the well known Schrodinger radial equation (SRE): 



2/i 



where \E'(r, 6', (/ 



[E - V{t)] - 



ui{r) = 



(6) 



Rni(r) Y, 



Imi y^i • 



is the complete wave function, r stands for the relative 



radial coordinate, and ui{r) = r Rni{r) is the reduced radial wave function. 

With respect to the boundary conditions, a regular solution near the origin for ui{r) 
could be |9] 

u{r 0) ^ r'+^ (7) 
Since asymptotically u{r 00) ^ we can take: 



u{r — >■ 00) — > exp 



2fi\E\ 



h 



(8) 



where \E\, as later will be seen, is an educated guess about the energy eigenvalue. 



The normalization condition reads 



dr\ui{r)[ 



dr\Ri(r)\ r 



2„2 



(9) 
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3. Numerical solution of SRE 



The SRE can be written as 

d?u{r) 



+ k{r-) u{r) = s{r) 



(10) 



Here (setting / = 0), k{r) = ^ [E — V{r)] is the kernel of the equation, and s{r) = 0. 



We can integrate these equations by means of the Numerov Algorithm [10] as follows: 

First we split the r range into N points according to r„ = r„_i + /i (where h is the step); 
then we write the wave function m„ = u{rn) = u{rn-i + /i), and A;„ = fc(r„) = fc(r„_i + h). 

Expanding u{r-) around r„: 



^3 ^4 

Un+i = u{rn + h) = u(r„) + W(r„) + yu"(r„) + yM"'(r„) + —u^'^\rn) + 0{h^) (11) 

h'^ 

Ur,-i ^ u{rn -h)= M(r„) - hu'{rn) + —u"ir^) - —u"'{rn) + TTrW^^H^n) + 0{h^) (12) 

Then approximating the second derivative by the three-point difference formula, and 
using it within the second-order differential equation we get the following recursive for- 
mulas, with a local error 0{h^): 

a) Forward recursive relation 



Ur 



2(1 - ^kn-l) - (1 + §fc«-2) Mn-2 

(1 + fj^n) 



(13) 



b) Backward recursive relation 

_ 2(1 - ^kn) Un - (1 + §fe«+l) Un+1 
(1 + f^^n-l) 



(14) 



Therefore, when we calculate our wave function using the backward-forward technique, 
we should note that the recursive formulas imply having knowledge of two initial values 
for each direction. 

It is also necessary to know the first derivative at the appropriate order. Following the 
above expansions, we then get: 



1 

2h 



:i + ^k. 



-l) Un+l - (1 + —kn~l) Un-l 




0(/l^ 



(15) 
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4. Wave function and normalization 

For the sake of simplicity, let us first focus on / = states. The SRE reduces to 

^ + p [i? - Vir)] u{r) = (16) 

To illustrate this point we choose an harmonic oscillator potential V{r) = /Jr"^, with j3 a 
constant. We then start, for instance, by using a forward calculation (with an appropriate 
energy eigenvalue). 

Since we are dealing with bound states, we find eigenf unctions at the classically al- 
lowed region with E > V{r) and classically forbidden region where E < V{r), they are 
separated at a turning point, r^, which can be estimated from the equality E = V{rc). 
If we perform a forward calculation, its asymptotic solution at the forbidden region may 
behave either as ~ e^"'" , where the positive value is non-physical. Thus, we have an 
admixture of those solutions and then, with successive iterations, the integration would 
be numerically unstable due to the dominance of the exponentially growing solution. As 
a general rule (TU], integration into a classical forbidden region tends to be inaccurate. 

Hence, for a given energy eigenvalue, we consider a calculation using both forward 
and backward solutions: from the allowed towards the forbidden region, with Uoutir) 
(outwards) eigenfunctions and from the forbidden towards the allowed region with Uin{r) 
(inwards) eigenfunctions . 

Let us note here that, to avoid numerical overfiows in the forward calculation, we do 
not usually start with u{r = 0): once included the centrifugal barrier term, the piece 
would originate an overfiow at r = 0. 



4.1. Bound state energy 



Since both Uoutij) and Uinir) satisfy an homogeneous equation, their normalization can 
always be chosen so that they are set to be equal at the Tc point. An energy eigenvalue is 
then signaled by the equality of derivatives at this point ^10]. At the matching point the 
eigenfunctions Uoutir) and Uin{r) and first derivatives u'^uti''^) Kni''^) must all satisfy 
the continuity conditions: 



Uout 



Ui 



U 



out 



(17) 



thus, we can write the corresponding condition for the logarithmic derivative at as 

(18) 



'^out 








- Uout - 


r-c 







and then we can define a G{E) function at whose zeros correspond to the energy 
eigenvalues as 



G{E) 



'^out 








-'^out- 


re 


-'^in- 





(19) 
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Therefore we proceed numerically in the following way: we set a trial energy range 
splitting this E range into N points, according to En = En^i + A^;, where A^; is the 
energy step. For each E^ we calculate their eigenfunctions Uout and at the Tc point; 
and we build the G{E) function here, looking for a change of sign in it (which implies a 
zero cross). Once we find it, we perform a fine tuning closing the energy range until the 
required tolerance. 



4.2. Matching eigenfunctions at the Tc point 

When we find the energy eigenvalue, the calculated inwards and outwards eigenfunc- 
tions will tend not to match at the point. However we can look for a strategy to solve 
this problem: 

Denoting the outwards and inwards functions directly obtained from the recursive for- 
mulas as $(r) and /(r), respectively, the physical Uguti''^) and Uin{r) eigenfunctions can 
be rewritten as 

Uoutir) = A$(r) Uin{r) = BI{r) (20) 

A and B are constants. Their respective derivatives are 
uUr)^A^'{r) u[^{t)^BI'{t) (21) 

By substituting eqs. (20) and (21) into eq.(17): 

and performing the difference, we get 



A 



B = feB (23) 



where /c will be a scaling factor to be applied to Uoutif). Therefore 
Uoutir) ^B *(r) u,„(r) - B I{r) (24) 

and S is a global factor that must be taken into account in the normalization process. 
4.3. Normalization 

Once the energy eigenvalue has been determined, we first insert it into the kernel, 
k(r), thereby generating their corresponding $(r) and /(r) functions; subsequently we 
calculate the factor. To find the remaining B factor, and therefore find the Ugut and 
Uin eigenfunctions, we use the normalization condition 

/ dr\ui{r)\^^l (25) 
Jo 
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where, following the asymptotic requirement u{r oo) 0, taking cutoff value. 

By separating the Uoutij) and Uinir) domains in the above integral, we can write 



dr\ui{r)\^ = I dr\uout{r)? + 



dr\uin{r){' 



(26) 



Using the equations (24) 



Jo 



rfr|M;(r)|2 = {Bf^f / rfr|<l>(r)|2 + B 

Jo Jr. 

the normalization condition then reads 



dr\I(r) 



B' 



{f,f drmr)\''+ rfr|/(r)|2 

Jo Jrc 



(27) 



(28) 



Denoting the result of the above integrals within brackets as A^, we can write B = 
thereby deriving the normalized eigenfunctions 



Uout{r) 
Uinir) = 



fc $(r) 



Kr] 



(29) 
(30) 



5. Integration and expectation values 



When performing the integration with the QQ-onia package, we use the following pro- 
cedure: If we name 



rN 



drf{r) ■ fn = f{rn) 



from the Euler-McLaurin summation formula [TT] with a given step h 



-N 



f + + + + + f 



B'2h'^ I I B<?,kh'^^ 
'{In ~ Jo) ~ 



2! 



(2A:)! 



where B2k are Bernouilli numbers. The first term of the r.h.s. in the above equality 
corresponds to the extended trapezoidal rule. 



If we set a number of steps to a multiple of 4, and apply the above formula for steps 
h, 2h, and 4h we obtain 



^ + /l + /2 + /3 + ... + f 



-2h 



2h 



y + /2 + /4 + ...+/7V-2 + ^ 



^'^^\fN-fo)-^^ifN-fo) + 0{h^) 



12 
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-4/i 



4h 



Y + /4 + /8 + ... + /iV-4 + Y 



12 4!30 

where T/v = Th = T2h = T^h- If we solve this system to ehminate derivatives up to 
0{hP) we arrive at the final formula of 

ft ^ ^ 64T, - 2OT2, + T4fe , ^,^6, 
c?r /(r) ^ — + 0[h ) (31) 

that is used in our calculations. 

Several expectation values are needed to be computed in our method. Once the eigen- 
functions have been normalized, one can calculate the expectation value of a given oper- 
ator O according to the definition 

f^max 

{0)= drul{r) O ui{r) (32) 

by using the above 0{h^) integration. O can be (if we want to obtain the mean square 



radius of the state y {r"^)), the potential or the derivative of the potential, etc. 
5.1. Square of the radial wave function at the origin 

We need to distinguish between the calculation of the square radial wave function at 
the origin (WFO) for Z = states, |i?„(0)p, and the calculation of the squared derivatives 
of the radial wave function at the origin for / ^ states, \R^^\Q)\^, with / = 1,2,3 
respectively corresponding to the P, D and F states. 

5.1.1. / = states 

From the well known calculation [3] derived from the Schrodinger equation we obtain 

\'^nimM? = ^ {V\r)) (33) 

where {V'{r)) is the expectation value of the derivative of the potential, and /i stands 
for the reduced mass. If we are dealing with hb or cc systems: yU = mQ/2 where mq is the 
heavy quark mass. 

For I = states the wave function is: "^nooii^) = Rno{i^), then 
\Rnm' = mQ{V'{r)) (34) 

Therefore, since the corresponding eigenfunctions have been already obtained, the eval- 
uation of the radial WFO reduces to a calculation that gives the expectation value of the 
derivative of the potential {V'{r)). 

Another method that can be used to evaluate the WFO is to extrapolate directly 
from the normalized eigenfunctions, taking values of the square |i?„(r)p = [ \u{r)\'^/r'^ ] 
from the region near to r = using an appropriate interpolating function. Even so, the 
previously stated technique yields a more accurate result, if the {V'{r)) calculation is 
reliable. Nevertheless both procedures can be employed together as a check. 
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5.1.2. I ^ states 

Here we consider the SRE, eq.(6), containing the centrifugal barrier term. If we rename 
the potentiah 



W{r) = V{r) + 



hH{l + 1) 
2/xr2 



the kernel of the SRE (see eq.(lO)) now becomes: k{r) = '^[E- W{r)]. 



1=1 CASE: The main goal here is to evaluate the square first derivative at r = 0, 
|i?^(0) p. Since we have a tool for directly calculating the derivatives u'{r), according 
to eq. (15), once we have the normalized eigenfunctions u{r), then using 



u{r) 



u'{r) — R{r) 



(35) 



we can extrapolate the square of the derivative, taking those values from the region 
near to r = with an appropriate interpolation function. 

1=2 CASE: We will obtain the second derivative at r = 0, |i?^(0)p but only in terms 
of its radial wave function and first derivative. To do this we proceed as follows: 
from SRE 



u {r) — — k{r) u{r) 



(36) 



If we calculate the second derivative using the above identity 



u{r) 



= - k{r) R{r) 



2 R'{r) 



(37) 



we can again extrapolate the |i?^(0)|^ function. The main advantage of this pro- 
cedure is that there is no need for any additional derivatives, and since u{r) and 
u'ij-) have already been calculated at the suitable order, the result holds with the 
foreseen precision. 

1=3 CASE: In the same way, using eq. (36), we can obtain the third derivative to 
extrapolate 



Rniir) 



uir] 



2 k{r) 



- k'{r) ) R{r) + ( ^ - k{r) ) i?'(r) 



(38) 



where k'{r) stands for the derivative of the kernel, k'{r) — 
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5.2. Calculating the heavy-quark velocity 

To obtain the ( f ^ ) value we perform a double calculation as a check: the first is from 
the Hamiltonian definition and the second uses the virial theorem. 

i) If r = f\—f2 is the relative radial coordinate between the quarks 1 and 2, with velocities 
\v\\= = its relative velocity v at the center of mass frame is 
V — 1v\ — —2v2- Then, we can obtain the quark velocity using the Hamiltonian, 
E — (T) + (y(r)), (where T represents the relative kinetic energy T — 




1 r 



E-{V{r)) 



(39) 



2fx I 



ii) Taking relative spherical coordinates, the virial theorem implies 



{T) = -{rV'{r)) 



Then, we derive the quark velocity from the expectation value of the product rV' (r) 
according to: 




1 



( r V'{r) ) 



(40) 



Apt 



in both calculations for bb or cc systems where 2/x = mg. 
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Table 3 



A J and Bj 


coefficients. 




j value 


Aj 


Bj 




-(16/3) 


+ 1 


'Pi 


-(4/3) 


+ (1/2) 


'P2 


+ (28/15) 


-(1/2) 



6. Bottomonia mass level splittings 



In this section we address different splittings of the mass levels of bottomonia in accor- 
dance with the expressions shown in section 2.1: 

n^Si — n^So splitting: In this case the only term that gives a non vanishing contribu- 
tion is the spin-spin term. Since from QQ-onia package we can previously calculate 
the nS state WFO, instead of eq.(4) we choose to use its final and well-known form 



A 



Mi'S^) - MCSo) 



a^(/i^ 



|i?no(0)P 



(41) 



where we will evolve the asiQ"^) to the appropriate scale [13 



where we know the quark-velocity (f ^) for each nS state from 



-oma. 



n^Pj and n^Pi splitting: We have to take into account altogether the spin-orbit and 
tensor terms, i.e. eqs. (3) and (5). Using a funnel potential, the explicit expressions 
to be computed read 



A 



M{'Pj) - M{centroid) 



Aj a^K" {^) + Bj ah'' {-) 



(42) 



where Aj and Bj are the corresponding coefficients for each case, as can be seen 
in Table 3; denotes again the evolved value up to the quarkonium scale. The 
expectation values are calculated using the wave function Vr corresponding to the 
nP centroid from a previous calculation, as we will see later. Once the ^Pj masses 
are calculated, the ^Pi value will be obtained according to [12], as 



M(iPi) 



5 M(3P2) + 3 M{^Pi) + M(3Po) 



(43) 
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7. Electric Dipole Transitions 

QQ-onia package allows us to calculate the wave functions of different states. Therefore 
one can calculate, e.g., the El transition rate initial {i) — ^ final (/) + 7 through the 
well-known expression 



TE^{^ --f + 1) = ^^^^ + ^ ' " ' ' ^^^^ 

where a = 1/137, and (eg) is the mean charge [ (—1/3) in the bottomonia case ]. In 
the QQ-onia file-example we will focus on ^5*1 —>-^ Pj radiative transitions with 5*,/ = 1. 
The photon energy, k is directly calculated from energy-momentum conservation law. If 
nif and mj stands for the experimental [H] masses of the final resonance (nPj) and the 
initial one (nS), respectively, one has 

2 2 

k = i (45) 

2 nii 

( / I r I i ) (in GeV^^) is the matrix element connecting final and initial state; we will 
evaluate it according to 



( / I r I i ) = j dr u*j{r) r Ui{r) (46) 

extracting the reduced wave functions of the nP and nS states from a previous calcu- 
lation with QQ-onia. 



8. The QQ-onia package 

The QQ-onia package is written with PAW software {Physics Analysis Workstation), 
which can be obtained for free from the CERN web site fTB] for several operative systems. 
This software contains a FORTRAN interface, called SIGMA. The QQ-onia package pro- 
vides a version of this software named PAW-NT for WINDOWS. The package also runs 
with UNIX-LINUX, pasting the bbnia-nl.kum files directly into any UNIX text editor. 

The package contains the files prepared to work with each of the previously mentioned 
cases. To illustrate how the machinery runs, and as a potential reference, we choose a 
standard static potential for heavy quarkonium: the Coulomb plus linear potential [6], 
and also as a further reference, we take a known set of parameters from |6], [16j with which 
to perform spectroscopy (for a comparison of the Cornell model with other approaches of 
the heavy quarkonium potential in the static limit, see |17j). 

The files solve the SRE for: / = states T{nS){n = 1, 2, 3, 4), and their file names are 
bbnia-ns.kum; the / = 1, x(nP)(n = 1,2) states, (bbnia-np.kum); and the / = 2, T{nD) 
ones, (bbnia-ld.kum). The package also contains a file called bbnia-4f.kum which is pre- 
pared for working with theoretical / = 3 [uF) states. The differences between these files 
is related to the centrifugal barrier and the calculation of either square of the radial wave 
functions or square of the derivatives at the origin. Since there are parts that are common 
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to both we will begin by explaining them. 

It must be pointed out here that in some parts of the software it is necessary to type in 
some analytic expressions (such as the derivative of the potential and the product rV'{r)). 
Although this procedure can be shortened, this is a good way to ensure control over the 
different variables, useful to analyze partial results. 

9. Getting ready to calculate energy levels: general settings 

First we open the file section SETTINGS in which we can read (to insert) the character- 
istic items: 

v/cr igina(l) r 0.9255 creates and defines the string tension (in GeV/ fm). 

v/cr mq(l) r 5. 18 creates and defines the b-quark mass. 

v/cr alf(l) r 0.39 creates and defines the strong coupling constant. 

v/cr hb(l) r 0. 19732858 is the he constant (in GeV fm units). 

ni sets the number of total steps; h=0.001 {fm) is the step. 

xO sets the minimum r value, rmm- the h value is the default. 

xin= [xO] + ( [ni] -1) * [h] sets the maximum r value, r^ax 

xc sets the cut r value Tcut according to Best = V{rcut), where we use Best as an estimated 
binding energy Best = Mexp — Sm^, and Mexp stands for the experimental mass of 
the resonance. Sometimes, if needed, we can also set rcut as the point at which the 
inwards integration has its first maximum. 

v/cr ele(l) r ' 'value' ' (with ele= 1,2,3) only appears in Z 7^ files, where it cre- 
ates and defines the eigenvalue of the angular momentum that needs to be inserted 
into the centrifugal barrier term. 

File section SIGMA APPLICATION BLCKl calculates the r range and defines the poten- 
tial: 

x=array ( [ni] , xO#xm) it calls si gma application (appl sigma) to establish the r range. 
The numerical values of ni , xO , xm, must be inserted here. 

POTENTIAL DEFINITION for / = files The V{r) value, (vo variable) is calculated for the 
full r-range. In our case vo=(igma*x)-(kfacl*(x**(-l))). 

POTENTIAL DEFINITION for / 7^ files ; 

vo=(igma*x)-(kfacl*(x**(-l)))+(bfac2*(x**(-2))), which contains the cen- 
trifugal barrier term, sets the W{r) value. 
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ENERGY SETTINGS The programme calls the sigma application (appl sigma) to es- 
tablish the energy range (using the variable e); here the numerical values of 
[ne] , E-min , E-mcLX must be inserted as follows: 

e=array (number of energy steps, minimum energy value# maximum energy value). 

ENERGY BLCK2 to insert the number of energy steps (ne). The programme will calculate 
the corresponding energy values, e, to cover the full energy range. 

BOUNDARY CONDITIONS this section of the programme runs automatically, the default val- 
ues arc tliosc previously discussed in section 2.2 u{r — > 0) ~ r'"*"^ and u{r — > oo) ~ 

j2ij.\B^,t\ 

exp —-^ — ^ r . 

those conditions are inserted into the two first values of Uouti''^) (uo variable) and 
Uin{r) (ui variable), respectively 

In summary ,we set h, [ni] , xO, xm, xc, [ne] , E-min, E-max;, the potential, and 
the corresponding constants. 

9.1. Steirting and determination of the bound state energy level 

Once we have inserted the previous settings, we write exec bbnia-nl .kum in the 
PAW interface. The programme then automatically searches for bound states within 
the [Emim Emax] range. If there is any change of sign in the G{E) function the pro- 
gramme stops, showing the new energy range [E'^^^, E'^^^^] in which the bound state can 
be found, and it asks for the next instruction 

(Type <CR> to continue or Q to quit). We then proceed in an iterative way: to im- 
prove the precision of the energy level, we skip the programme (by typing in q) and on 
opening the file, we type in the new energy range in the ENERGY SETTINGS file section. 
We then restart the programme and repeat this process until the energy level has the 
desired tolerance. 

Once the programme stops showing a suitable energy range [El^^l, El^^l], we push the 
enter key. The programme will then ask for the final (ee(l)) energy eigenvalue, then we 
type in this value (which could be the mean value of -B^a* ]), and it will calculates 

automatically the eigenfunctions Vr, showing at the screen the normalization proof. 

As a final result we get: 

1 The normalized reduced radial wave function, u{r) throughout the full range, stored as 

a [ni] -dimensional vector (yl) 

2 The u*{r) u{r) product {in fm~^ units)'which is also stored as a [ni] -dimensional 

vector (f c). These values can be either displayed in a graph, by typing in v/dr f c, 
or obtained numerically (v/wr fc). 
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9.2. Mean squEire radius 

File section **BEGINS SqR(<r2>) CALCULATIDN** . The software automatically calcu- 



lates the mean square radius of the state {y {r"^)) according to its expectation value. The 
[ni] -dimensional variable (f i) to integrate is f i=fc*xx, where xx is the variable. As 
this is the final result, this value is stored at the rad variable (in fm units), and displayed. 

9.3. Computing the heavy-quEirk velocity using the virial theorem 
File section <v2>/c2 calculation USING VIRIAL THEOREM** 

1 The first step is to calculate the expectation value of the product of X- W, ( ( rV'{r) ) ); 

to do this, BEFORE running the programme, we must TYPE in its corresponding 
expression at the variable xdpot. Using our example we must insert: 
xdpot=(igina*x) + (kf acl* (x** (-1) ) ) 

2 The programme integrates the variable f i=f c*xdpot and performs the velocity calcu- 

lation. The result is displayed and stored at the variable v2. 

9.4. Computing the heavy-quark velocity using the Hamiltonian 
File section **2) <v2>/c2 calculation USING <H>=<T>+<V(r)>** . 

1 The programme first calculates the expectation value of the potential {V{r)). The 

variable to integrate is ff i, where ff i=vo*f c. The result is stored at the variable 
vbar. If desired, it can be read by dropping the comment variable * at line: 
*v/pr vbar. 

2 It calculates the velocity value using the variables: vbar, the energy eigenvalue, ee, 

and the inverse of the quark mass, inq, and displays the final result. 

10. Computing the radial wave function squeired and derivatives at the origin 

Squeire of the radial WFO (For I — states: bbnia-ns . kum files). 

File section **BEGINS WAVE FUNCTION AT THE ORIGIN CALCULATION** . First the 
expectation value of the derivative of the potential {V'{r)) is calculated; to do this, 
BEFORE the programme runs, we must type in its corresponding expression at the 
variable dpot. Using our example, we must insert: 
dpot= (igma) + (kf acl* (x** (-2) ) ) 

The variable to integrate is f i, where f i=dpot*f c. Later, the final calculation is 
performed and the result is displayed (in GeV^ units) and stored at the wf o variable. 

SquEire of the first derivative of the radial WFO (For I — 1 states: 
bbnia-np .kum files). 

File section WAVE FUNCTION DERIVATIVE CALCULATION . General settings are equal 
than before by changing V{r) — >■ W{r), then the centrifugal barrier term appears 
(bfac2*(x**(-2))). As a result we obtain, according to eq.(35), the squared first 
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derivative of the radial wave function |i?^(r)p for all r value (in GeV^ units). The 
result is stored at the [ni] -dimensional variable derc. Then, to obtain its value 
at the origin, we must export the (derc) values in order to to extrapolate with an 
appropriate tool (such as the PAW inner routine vector/fit). By typing in v/dr 
derc, we can plot these values to select the range before exporting. To do this 
numerically type in v/wr derc. 

Squaire of the second derivative of the radial WFO (For I — 2 states: 
bbnia-nd.kum files). 

File Section WAVE FUNCTION DERIVATIVE (2) CALCULATIDN We obtain, by using the 
eq.(37), the second derivative of the radial wave function (in GeV'' units) |i?^(r)p 
values for all r, these values are stored at the [ni] -dimensional variable sed2, and 
ready to be extrapolated. 

Square of the third derivative of the radial WFO (For theoretical I — 3 states: 
bbnia-nf .kum files) 

File section WAVE FUNCTION DERIVATIVE (3) CALCULATION . Before starting the pro- 
gramme, we must TYPE in the corresponding derivative W'{r) expression into the 
variable kerdO. Then the programme will calculate the derivative of the kernel, 

k'{r). As a result, according to eq.(38), we have the [ni] -dimensional variable 
thrd2, which stores the full r-range of the third derivative of the radial wave func- 
tion K'(r)|2 (in GeV^ units). 

10.1. Wave function plot 

When a bbnia-nl file stops, a plot of the dimensionless reduced wave function can 
be obtained (from the same PAW screen), by running the routine graph-nl.kumac. We 
obtain the plot of the u2 variable as output, which is the [^(r)^ value multiplied by the 
Bohr radius of the resonance, Oq = ^/{Cf //). 

11. QQ-onia package spin-dependent and El 

This part of the software is organized within QQ-onia package as follows: The file 
SSplit-nS.kumac allows to calculate the n^Si — u^Sq splitting for each n level. The files 
Split-nP .kumac {n = 1,2) calculate as example the n^Pj and n^Pi splitting. There is 
also a file E1-2S IP. kumac, which analyzes the El 2S — > IP transitions, which has two 
subroutines named ini-2s.kum and final-IP. kum devoted to generate, respectively, the 
initial and final states for the matrix element calculation. 

SSplit-nS.kumac This file has mainly three blocks: the first one ask for the values of 
the quark mass (the mq variable), the value of the wave function at the origin (wf o 
variable), and the quark velocity of the nS state (v2q variable), whose values can 
be previously found with the bbnia-ns.kum files. Once these values are entered the 
programme runs automatically; it calculates the quarkonium scale and the ag value 
(alf variable) through the alpha-s Evolution Block. Later using the cq.(41) 
(the Delta-ss=[M(n3Sl)-M(nlS0)] calculation block) it displays the result of 
the energy splitting (in GeV) through the variable deltss. 
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Split-nP.kumac The first part of ttie file generates the centroid wave function, (the 
GENERATION OF THE CENTROID WAVE FUNCTION-BLOCK. This implies to run first 
the corresponding bbnia-np.kum in order to know the energy eigenvalue, the r^ax 
and Tc values to be inserted here [We insert the energy eigenvalue (ee variable) 
when the programme ask for it (EIGENVALUE? (GeV))]. This block can be performed 
in other cases by pasting from the bbnia file the part corresponding to the block 
labeled ONCE ENERGY LEVEL FOR BOUND STATE IS DETERMINED*************** 
WAVEFUNCTION CALCULATION. At the end of this part we have stored the normalized 
wave function within the yl variable. Later the programme will run automatically, 
it calculates the < > and < > expectation values, evolves the as, and 
applies the final expression eq.(42). Finally it displays the energy differences and 
the final ^Pj masses [Also the ^Pi mass using eq.(43)] by means of 

M(^Pj) = AEj + Mcentroid 

where AEj is the calculated energy difference of each J and Mcentroid is the 
centroid mass previously calculated from the bbnia-np.kum file. 

El-2SlP.kumac This file, as a complementary tool, is an example of how to handle 
QQ-onia with the El transition rates between 2S and the ^Pj states. First, the 
programme calls the ini-2s.kum subroutine in order to create the initial 25" state, 
which is performed by inserting the results found using the bbnia-2S.kum file: the 
energy eigenvalue (at the Enter the energy eigenvalue line), the Vmax and 
values; then it generates its corresponding wave function. Later, the software calls 
the final- Ip.kum subroutine and, in the same way, it creates the wave function 
corresponding to the final state. It must be pointed out that, when we want to 
modify the common parameters used generating the final and initial states, we 
perform it within the El-2SlP.kumac file. The programme then runs automatically: 
it calculates, using eq.(46), the matrix element < f\r\i > (stored at the me variable 
in GeV~^), the photon energy for each final J (the kph0,l,2 variables, in GeV) 
according to eq.(45), later it computes the final values, in keV, by means of the 
eq.(44), the result is displayed for each J = 0, 1,2 value through their respective 
gmaO ,1,2 variables. The remaining factors and settings are explained and displayed 
inside the El-2SlP.kumac file. 

12. Numerical accuracy 

One important issue is determination of the numerical accuracy of the results. Since 
no complete bottomonium wave function has yet been available, we can make our check 
using an harmonic oscillator potential, which has a known analytical solution, but it is 
numerically sensitive (as previously mentioned, this potential exhibits numerical instabil- 
ities due to their exponentially growing solutions). It is therefore a good candidate for 
checking our software. 
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We test V{r) = \ix uj'^ by inserting it in the file block SIGMA APPLICATION BLCKl, 
and by inserting in the SETTINGS file section /i = 100 MeV/c^; Tic = 197.32858 MeV/m; 
u = 2 fm~^c. As boundary conditions ^ we use u{r — 0) ~ r'"*"^ and for u{r — >■ oo) ~ 
e~^^ , where /? = 

The analytical energy levels are E{n,l) = (2n + l + l^huj. 

The reduced wave functions expressed in terms of the generalized Laguerre polynomials 
[9] are: 



Uni{r) 



T{n + 1 + 3/2) 



1/2 



(47) 



We then analyze levels IS* = (0,0); 25 = (1,0) and 2D = (1,2). The results can be 
found in Table 4. 



Table 4 

Harmonic oscillator results. 



Energy level (MeV) 


Analytic 


Numeric 


E(0,0) 


591.986 


591.999 


E(1,0) 


1381.330 


1381.315 


E(2,l) 


2170.614 


2170.628 



(0,0) LEVEL Taking 

Tcut as the turning point (using £"(0,0) = Vircut))-, we find that 
"^cut = 1-72 fra with Vmax = 5 fm and r^m = ^, with a step h = 0.01 fm, thus 
[m] = 500. We start searching for the associated energy level through the range 
E{Q, 0) e [400, 800]MeV, with an energy step of A^j = 1 MeV, i.e. [ne] = 401. 

(1,0) LEVEL We repeat the procedure, but now Vcut = 2.6 fm with Vmax = 5.6 fm, thus 
[ni] = 560; starting with an energy range E{1, 0) G [1100, 1500]MeV, and an energy 
step oi Ae = 1 MeV. 

(1,2) LEVEL First the centrifugal barrier term ^-^^ must be added to SIGMA APPLICATION 

BLCKl. We use Vcut = 1-08 fm, r^ax = 6 fm, with an energy range 
E{1,2) E [1900,2300]MeV, A^^ = 1 MeV. 

Concerning the reduced radial wave functions u{r)ni, in all cases, we find a deviation of 
less than 1% with respect to the exact value from eq.(47). Figure 1 shows the result for 
the u{r)i2 case; Figure 2 shows its corresponding A|u(r)|^2/I'"('")li2 relative error (in %). 
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Figure 1. Harmonic oscillator (n = 1,/ = 2) reduced wave function. 



13. An example: Numerical Results 

To illustrate how to the QQ-onia package runs, all attached bbnia-nl.kum files can be 
opened; their corresponding settings are listed inside. Here, we will focus on the T(IS') 
case (bbnia-ls .kum file) in order to explain some relevant details. 

We look for the lowest state of the bb family. First we estimate its energy eigenvalue, 
Bl^t, from the experimental data [H], then B^^^ = M^^^ - 2mb = -0.8997 GeV, (with 
^exp ~ 9.4603 GeV and mj, = 5.18 GeV). To ensure that we find the lowest state we set a 
wide energy range below the B^^^, thus we start typing at ENERGY SETTINGS section a trial 
energy range B^^ G [— 5.,0.] GeV with an energy step of A^; = 0.5 GeV, i.e. [ne]= 11. 

Using = V{rcut), we find Vc ~ 0.1 fm (alternatively, instead of Best = V{rcut), we 
can take the Tc point where the inwards integration has its first maximum). To ensure 
the asymptotic conditions we then take Vmax = 1-0 fm; taking h = 0.001 fm, the number 
of steps is [ni]= 1000. 

We type in exec bbnia-ls . kum at PAW screen, then the programme runs until it stops 
displaying: 

bound state around 
E(10) = -0.5 
E(ll) = 

Type <CR> to continue or Q to quit 

This is the new energy range (in GeV) within which we can find the eigenvalue of the 
bound state. We then stop the programme (q) and insert these new values at the ENERGY 
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(iuVu') (percent) (harmonic oscillalor(l ,2)) 



100 200 ^00 400 



Figure 2. Harmonic oscillator (n — 1,1 — 2) relative error. 



SETTINGS section and restart the programme. We repeat the procedure, and after five 
iterations the screen displays: 

bound state around 
E(10) = -0.170209 
E(ll) = -0.170208 

Type <CR> to continue or Q to quit 

As there is sufficient precision, when we press the enter key, the screen shows: 

BOUND STATE ENERGY? 
EE(1) 

We then take the mean value of the above quantities —0.1702085 (GeV) as the eigen- 
value, and type it in at the screen. 

EE(1) -0.1702085 

As the calculated eigenvalue Bl^^^ docs not match the Blf^ value, we can establish a 
scale factor (F) to be applied to the whole spectrum: 

F = Mll^ - (2m, + S^^i) 

(Alternatively, we can redefine the quark mass according to 2m'j, = 2m;, + F) . The 
masses of the higher states M^'^^ can therefore be obtained from 

M^L = ^rn, + B:Ii, + F 
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where, -BcaZc stands for the calculated energy eigenvalue obtained from the software for 
each higher state. 

In Table 5 the results for the masses of several resonances can be found. 
Table 5 

Coulomb plus linear potential: mass results (in GeV units). 



bb LEVEL Experimental mass Mass from [6] Mass from QQ-onia 



T{IS) 


9.4603 


9.4603" 


9.4603" 


x{iPY 


9.9001 


9.96 


9.9584 


T{2S) 


10.02326 


10.05 


10.02772 


T{ID) 


10.1622 


10.20 


10.2080 


x{2Pf 


10.2620 


10.31 


10.3125 


T(35) 


10.3553 


10.40 


10.3971 


bb{AFY 






10.3995 


T{AS) 


10.5794 


10.67 


10.6739 



(a) used to set the ground level in both references 
(6) Xbj{nP) centroid (c) Theoretical level 



After inserting the energy eigenvalue EE(1), the programme will calculate the remain- 
ing quantities automatically: it first shows the normalization check and then displays a 
plot of the f c value. 

normalization proof 
NNOR(l) = 1 

It continues to show the value of the square of the wave function at the origin 
RADIAL WAVEFUNCTION AT THE ORIGIN (in GeV^) 
WFO(l) = 14.0927 

If we run a / 7^ file, we do not see this screen but when the programme finishes we 
can export the corresponding \R^\r)\'^(yr) and extrapolate. 

Table 6 summarizes the results for the square of the WFO (or its derivatives). 

After the WFO value, the calculation for the mean square radius calculation from the 
IS* file will show 
SQR(<r2>) (in fm) 
RAD(l) = 0.201043 



The results for each resonance can be seen in Table 7. 
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Table 6 

Coulomb plus linear potential: |-R^'^(0)p values in (GeV)*^'^^^'^ units. 



hh LEVEL 


|i?(,')(0)p from QQ-oma 


\R\^>{0)\^ from ^ 


T(l^) 


14.09 


14.05 




2.062 


2.067 


T(25) 


5.947 


5.668 


T{ID) 


0.835 


0.860 


X(2P) 


2.440 


2.438 


T(35) 


4.276 


4.271 


hh{AFY 


0.551 


0.563 


T{AS) 


3.675 


3.663 



Table 7 



Coulomb plus linear potential: mean square radius (in fm). 



hh LEVEL 


\J ir"^) from Q 


Q-onia 


ir"^) from [B] 


T(15) 


0.20 




0.20 




0.38 




0.39 


T{2S) 


0.46 




0.48 


T{ID) 


0.52 




0.53 


X(2P) 


0.63 




0.64 


T(35) 


0.71 




0.72 


hh{AFY 


0.64 






T(45) 


0.91 




0.92 



Finally, the IS file displays the velocity result obtained applying the two methods 

(<f^>/c^) [virial theorem] 
V2(l) = 0.0962335 
(<u2>/c2) [<H>=<T>+<V>] 
V2A(1) = 0.0962476 

and the corresponding results for each case can be seen in Table 8. 

Finally, as an illustrative example, one can observe in Figure 3 the dimensionless T(IS') 
reduced wave function taken from the graph-nl.kumac file. 




Figure 3. T(IS') reduced wave function. 



14. More results using another potential 

As a further example, we present results obtained using a Leading Order potential for 
a heavy quarkonia system 

T/(l) 

nib 

which contains a V^^^ static term (which could be a Coulomb plus linear potential), and 
an additional 0{l/m) piece whose contribution is comparable to the static part. From 
Lattice analysis p^HTO] , we set the explicit form of this potential as follows (the so called 
Cornell-modified potential, see [20] and references therein for a complete description): 

c c' 

Vcor-mod{r) = -- - ^ + (TV + fl (48) 

The values of parameters m^, c, c', a and were obtained by applying a fitting procedure 
to the bottomonium spectrum (using T(IS') and T(2S') states), not from Lattice estimates. 
From the fit we obtained the following values for the parameters of the potential: 
TJib = 4.7 GeV, (/i = 0) 

a = 0.217 GeV^, c = 0.400, c = 0.010 GeV"^ 

Here, the quark mass redefinition according to 2m'{, = 2mb + F was used. 

The values of the predicted masses, square of the WFOs (or their derivatives) , and 
other parameters of interest for different bottomonium states obtained using this potential 
are shown in Table 9. Its it possible to observe an excellent agreement with experimental 
mass values. 
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Table 8 



Uouiomb plus 


liner potential: velocities { v"^ ). 




00 LhvVbL 


\ V ) ironi 


Q-onia"" 


( V ) irom [6j 


T(l^) 


0.096 




0.096 




0.065 




0.065 


T(25) 


0.078 




0.076 


T{ID) 


0.067 




0.067 


X(2P) 


0.075 




0.076 


T(35) 


0.085 




0.085 




0.073 






T(45) 


0.096 




0.097 



a : Results doubly checked from Hamiltonian and Virial methods. 



Table 9 

Values of the predicted and experimental masses (in GeV), square of the WFO 
(or derivative) in GeV^^^', mean square radius (in fm) and the typical quark velocity 
for the T(1S', 25", 35", 45'), Xb(l-P; 2P) and T(1D) states for the Cornell-modified potential. 



Resonance 


Mass 


Exp. 








T(15) 


9.4603 


9.4603 


12.65 


0.23 


0.094 


X5(1P) 


9.8929 


9.9001 


1.409 


0.40 


0.071 


T(25) 


10.0236 


10.0233 


6.444 


0.51 


0.091 


T(1D) 


10.1476 


10.1622 


0.562 


0.53 


0.078 




10.2729 


10.2600 


1.854 


0.63 


0.089 


T(35) 


10.3750 


10.3552 


5.404 


0.71 


0.103 


T(4S) 


10.6477 


10.5794 


5.194 


0.88 


0.120 



15. Spin-dependent splittings and E\ width: Numerical Results 

In this section we show some results from the above mentioned files. We would like to 
stress that such results are not exhaustive, but just to illustrate how QQ-onia runs. 

n^Si — n^So splitting: In Table 10 we summarize the results obtained for the n^Si — 
n^So splitting running the S Split, kumac file using different potentials. According to 
section 6, the parameters involved in this calculation are: the heavy-quark mass, 
its corresponding velocity, and the radial wave function at the origin (rWFO). We 
present our results from QQ-onia using first the set of parameters from the Cornell 
potential [6]. We also calculate the splitting employing the parameters from the 
Buchmiiller and Tye (B-T) QCD-motivated potential [21]. Lastly, we have reob- 
tained as a final check the hyperfine splitting of a Coulomb plus Linear Potential of 
Ref. [22j (for the particular case u = 1) . 
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Table 10 



AM = M{^Si) 


— M{^So) values 


(in MeV). 






nS level 


Cornell 


AMb-t 


AMcpL-type 


A Experimental 


IS 


135 


74 


66 


71.4 


2S 


57 


38 


44 




3S 


42 


28 


33 





Parameters employed [|i?„o(0)P(GeV'^) and mfe(GeV) ]: 
Cornell: m,, = 5.18, |i?io(0)P = 14.06, |i?2o(0)P = 5.668, |/?3o(0)P = 4.271 
B-T: rUb = 4.88, |i?io(0)P = 6.477, |i?2o(0)P = 3.324, |i?3o(0)P = 2.474 
CpL-type: = 5.1, |i?io(0)P = 6.173, , |i?2o(0)P = 4.027, |/?3o(0)P = 3.080 



A comment is in order here. Recently the BaBar collaboration has claimed the 
discovery of the long-awaited rjb{lS) state [23]. Its observed mass (~ 9.389 GeV) 
is somewhat lower than expected yielding a somewhat large hyperfine splitting. 
As can be seen in Table 10, the Cornell potential provides the largest deviation 
(by excess!), obviously due to the fact that this potential yields the largest value 
for the rWFO, while the rWFO values from the B-T and CpL-type potential are 
considerably smaller. In fact, let us note that the BuchmuUer-Tye model provides 
an excellent agreement with the experimental result. 

Table 11 



ni^Pj) and n{}Pi) Masses (in GeV); Cornell potential. 



nS level 


QQ-onia 


Experimental 


l^Po 


9.867 


9.860 


l3Pi 


9.893 


9.893 




9.900 






9.911 


9.913 


2'Po 


10.236 


10.232 


2^Pi 


10.256 


10.255 


2^Pi 


10.262 




2'P2 


10.271 


10.269 



n^Pj and n^Pi splitting: In Table 11 we summarize the results obtained from the 

Split-nP.kumac files using the Cornell potential, in accordance with eqs.(42, 43), 
to be compared with the experimental data from [H]. 

El 2(^Si) — > l(^-Pj) Transitions: To illustrate this point we show in Table 12 the result 
obtained from the El-2SlP.kumac file, in accordance with eq.(44). The Cornell 



29 



Potential was employed to generate the initial and final state wave functions. A 
nice agreement with experimental data is found. 

Table 12 

El [2{^Si) l^Pj)]. From QQ-oma < f \ r \ i >= 1.6915 GeV"^ 



Final state 


Photon energy (MeV) 


Width (keV) 


Experiment (keV) 


J = 


162.48 


1.47 


1.22 ±0.16 


J= 1 


129.63 


2.25 


2.21 ±0.22 


J = 2 


110.44 


2.32 


2.29 ±0.22 



Experimental data obtained from [H] . 



16. Summary 

The main goal of this work is to provide a multipurpose (user-friendly) package to ob- 
tain the wave functions at the origin and other relevant properties of heavy-quarkonium 
systems, assuming a basic knowledge of the PAW software by the interested reader. 

Besides, we would like to stress some special peculiarities of our package, such as pro- 
viding an easy procedure to normalize the resulting wave functions on account of the the 
Numerov forward-backward framework. In addition, the calculations of the /-derivatives 
at the origin for angular momentum I = 2,3, ... only require a numerical computation of 
the wave function and first derivative through the analytic expressions given in the main 
text, thereby keeping the suitable precision for higher derivatives . 

In addition, we present an alternative method to estimate the heavy-quark velocity 
using the well-known virial theorem and a quick 6*''-order integration, which can be con- 
sidered as well as a check of the traditional calculation. In fact, this method can also 
be interpreted as an indirect test of the goodness of the potential-probe together with 
the (whole range) prior calculated wave function, since it uses the expectation value of r 
times the derivative of the potential { r V'{r) ). 

Moreover, some worked examples are presented for the bottomonium system using a 
Cornell-type potential. Finally, an additional machinery has been implemented in the 
code containing files to analyze the impact of the spin-dependent terms in the potential, 
as well as a tool dealing with El transitions. Another set of worked examples is presented 
in this regard. 
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